Spin-resolved spectra of Shiba multiplets from Mn impurities in MgB2 
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We study the effect of magnetic Mn ions on the two-band superconductor MgB2 , and compute both 
the total and spin resolved scanning tunneling spectrum in the vicinity of the magnetic impurity. We 
show that when the internal structure of the Mn ion's d-shell is taken into account, multiple Shiba 
states appear in the spectrum. The presence of these multiplets could alter significantly the overall 
interpretation of local tunneling spectra for a wide range of superconducting hosts and magnetic 
impurities. 
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I. INTRODUCTION 

The interaction between a single magnetic impurity 
and the superconducting host reveals fundamental prop- 
erties of both the magnetic ion and the host material. 
This interaction was first studied theoretically, within 
the framework of BCS superconductivity. During the 
late sixties Shiba^ showed that a magnetic impurity pulls 
down from the continuum states a pair of bound states 
inside the superconducting gap. Indirect indication for 
the presence of finite spectral weight inside the gap of 
an impure superconductor could be inferred from global 
probes of the density of states. However, direct evidence 
for the existence of the so-called Shiba states requires 
an accurate measurement of the local density of states 
near the impurity. Such a measurement became avail- 
able only recently by using high vacuum, low temper- 
ature scanning tunneling spectroscopy (STS). Yazdani 
and his coworkers imaged^ the local density of states 
around Mn and Gd impurities deposited onto Nb single 
crystals. They found clear evidence for localized states 
in the vicinity of the magnetic impurities, in qualitative 
agreement with Shiba's original findings, and also with 
their own model calculation based on a non-selfconsistcnt 
solution to Bogoliubov-de Gennes equations. Quanti- 
tative discrepancies, however, are also clearly present, 
especially when comparing the width and spatial de- 
pendence of the resonances to theoretical expectations. 
The presence of magnetic impurity induced bound states 
in a superconductor was turned around and used, both 
theoretically- '"^'^ and experimentally^, as an investigative 
tool to probe the unusual ground state of the cuprate su- 
perconductors. 

Although there exist some precious numerical renor- 
malization group and Monte Carlo results for quantum 
dots attached to superconducting electrodes,^'* most of 
the theoretical studies carried out so far for magnetic im- 
purities in a superconductor follow Shiba's original work, 
and use predominantly a classical spin model to describe 



the magnetic impurity and assume a single spin one- 
half electron channel that couples to the magnetic impu- 
rity. Furthermore, the coupling is assumed to be in the 
s-wave channel, and spin-orbit coupling is generally ig- 
nored. This set of approximations worked beautifully for 
most of the experiments performed so far and provided 
simple, elegant and intuitive results. However, recent ad- 
vances in the resolution, stability and processing of scan- 
ning tunneling imaging opened the door for visualizing 
structures that go beyond the class of Shiba-like models. 
Indeed, magnetic impurities have a more complicated in- 
ternal structure^: The magnetic moments are usually due 
to low-lying and crystal-field split d- or /-levels with mul- 
tiple occupancy. The aim of this paper is to demonstrate 
that (i) the internal structure of the Mn impurity has a 
major impact on the structure of the Shiba states, and 
(ii) these novel features should be readily observable with 
the current resolution of STS measurements. In partic- 
ular, multiple channels of charge carriers couple to the 
magnetic impurity through channel-dependent coupling. 
The combination of these ingredients generally leads to 
the appearance of multiple pairs of Shiba states. We com- 
pute the spatial and spin structure of the scanning tun- 
neling microscopy (STM) spectra around the magnetic 
impurity and show that these states appear as distinct 
resonances inside the superconducting gap, and can be 
most clearly resolved in spin resolved STM spectra. 

In the following we illustrate our results on the spe- 
cific case of Mn-doped MgBj , but we wish to emphasize 
that much of our discussions carry over to other systems 
as wellAfi', and that our conclusions are rather general. 
There are several reasons to choose the Mn — MgBj sys- 
tem. Despite the relatively recent discovery of its es- 
sentially conventional superconducting phase, MgB2 has 
been thoroughly characterized both experimentally and 
theoretically^^, and therefore provides an ideal testing 
ground for our theoretical framework. Several materials 
parameters of MgB2 are also in a convenient range for our 
investigation. First, in order to observe a Shiba state by 
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scanning tunneling spectroscopy (STS), one needs a rel- 
atively large gap. MgB2 is a perfect candidate in this 
respect since it is a conventional superconductor that 
has an unusually high critical temperaturei^, — 39 K. 
Second, MgB2 has a hexagonal AlB2-type structure and 
a highly anisotropic band structurei^ii^. As we shall 
see below, this leads to a clear separation of the mul- 
tiple Shiba states. The presence of two gaps in MgB2 
has been well established by now through a variety of 
spectroscopic probe o^^i^^i^^'^^ . It is therefore an interest- 
ing question, how the presence of these two gaps influ- 
ences the structure of Shiba states. Although a series of 
experimenta l^^i^° and theoretica P^i^^ investigations have 
been recently completed for MgBg doped with nonmag- 
netic as well as magnetic impurities, no experimental or 
theoretical study has been completed for the local elec- 
tronic structure of a single magnetic impurity in this 
compound. This paper now provides a detailed theoret- 
ical discussion of the single magnetic impurity problem 
in MgB2 and other superconductors where the multiple 
degrees of freedom of the conduction electrons and the 
impurity could lead to experimentally observable conse- 
quences. 

Finally, there is another advantage for studying Shiba 
states in MgB2. The strong coupling, short coherence 
length and the consequently robust condensate allows us 
to investigate the effect of a single magnetic impurity 
on a superconductor in the regime where the order pa- 
rameter remains spatially constant. Here the results of 
Flatte and Byers** are quite valuable: According to their 
calculations, for a superconductor with coherence length 
^kp = 10 the relative spatial fluctuations in the local 
order parameter remain below 5%, even at the impurity 
site. This result is valid for the entire range of interest 
< 5 < 1 for the dimensionless coupling g (between the 
magnetic impurity and the superconducting quasiparti- 
cles, see below). As shown in the following sections of 
this paper, the spatially constant order parameter pro- 
vides considerable simplifications in our calculations, and 
this model allows us to make experimentally testable pre- 
dictions for the presence of the multiple Shiba states in 
MgB2. 



II. HAMILTONIAN 

A. Band structure calculation 

As mentioned above, MgB2 crystallizes in the hexag- 
onal AlB2-type structure^'^ in which the ions consti- 
tute graphite-like sheets in the form of honeycomb lat- 
tices separated by hexagonal layers of Mg ions. Band 
structure calculations^^ indicate that Mg is substantially 
ionized, and the bands at the Fermi level derive mainly 
from Boron p orbitals. Four of the six p bands cross the 
Fermi energy, and the Fermi surface consists of quasi-2Z? 
cylindrical sheets, due to B - p^.y orbitals, and a 31? tubu- 
lar network (mostly originating from B - orbitals). It 



is believed that both structures participate in the forma- 
tion of the superconducting state, though the gap is very 
different on the tubular network and on the cylindrical 
sheets. 

Let us first discuss the tight-binding Hamiltonian we 
use and the corresponding band structure. In spite of its 
simplicity, this tight binding description is rather robust, 
as can be checked by a direct comparison to the results of 
more sophisticated ab-initio band structure and density 
of states (DOS) calculations.^^ In the rest of the paper 
we shall use the following simple Hamiltonian to describe 
the normal state of MgB2 , 

(*i,a,a*r',a'.a + /i.e.) , (1) 

where /.t sets the Fermi energy and ^r.a,(7 is the annihi- 
lation operator of an electron of spin a on p-orbital a 
(a = Px,Py,Pz) of the B ion at position r, 

r = R + d . (2) 

The vector R in this expression points to the center of 
the unit cell and d gives the position of the B ion within 
the unit cell. Note that there are two atoms per unit 
cell, which shall be labeled by the index 6 = 1, 2 in what 
follows. The hopping matrix elements t"'^" in Eq. ([T|) 
connect only neighboring sites, but their value depends 
on the relative orientation of the p-orbitals. Quasiparticle 
energies are measured from the Fermi energy, fi. 

The Hamiltonian above can be easily diagonalized in 
Fourier space. The field operators '^T,a,a can be ex- 
panded as 

■^r.,a,a = ^ E ^ ^k,f,,a , (3) 

k,& 

where fl is the number of unit cells, and Ck,b,cr is the an- 
nihilation operator of an electron in band b (b — 1 , . . . , 6) 
with momentum k, spin a, and energy ek,b- The band 
energies and the wave function amplitudes eb-^a,S are de- 
termined by the eigenvalue equation 

Ha,S-a',S' (k) eb-a',S' — £k,b eb-a,5 , (4) 

q' ,5' 

where i?a,(5;a',(5'(k) is essentially the Fourier transform 
of the hopping matrix, detailed in Appendix A. In our 
tight binding model we have six bands: Four of them 
derive from p^y orbitals while the remaining two from 
Pz orbitals. The band structure obtained is presented in 
Fig. [1] Notice that both pz (tt bands) cross the Fermi 
surface but only two of the px,y bands (cr bands), cross 
it. 

In the presence of superconducting order, one must 
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FIG. 1; Band structure for MgB'z, along the symmetry lines, 
computed in the framework of tight-binding model described 
in Appendix A. There are six bands: four from px,y orbitals 
(solid lines) and two from (dashed lines). Both bands 
and only two px,y band cross the Fermi level which corre- 
sponds to zero energy, Ef = 0. 



modify the Hamiltonian above and add the pairing terms. 



fc,k 



Here the summation goes over those four bands that cross 
the Fermi energy (6 = 1 ... 4) . We assume further that 
the superconducting gaps take only two different values: 
in the Px,y bands Axy ~ 7.5meV while for the p^-bands 
it is Az « 2.5 meV. In our work, we shall neglect fur- 
thermore the position-dependence of the gaps around the 
magnetic impurity. This approximation is justified by the 
short coherence length in MgB2 , as already explained in 
the introduction. 



B. Interaction with a magnetic impurity 

To carry out a quantitative analysis of the magnetic 
impurity problem, we first need to establish how mag- 
netic spins couple to the conduction band. The interac- 
tion part of the Hamiltonian depends on the specific lo- 
cation and electronic structure of the magnetic impurity 
considered. In what follows, we provide a detailed anal- 
ysis for Mn impurities, which have already been doped 
into MgB2 , though similar considerations hold for other 
types and positions of magnetic impurities. Mn ions pre- 
sumably substitute the Mg atoms, and most likely take 
an Mn^^ configuration with a half-filled d-shell and a spin 
S « 5/2i^ As shown in Fig. [21 the five- fold degeneracy of 
the d-states is lifted by the local hexagonal crystal field 
into three multiplets that we can label by the original 



angular momentum quantum numbers /i of the d-states, 
1^). Each of these states is occupied by a single electron, 
and hybridizes through a hybridization with a specific 
local combination of p-states, V'/i that we construct next. 

The Mn ion is in the middle of a cage of 12 B ions, 
that we shall label by the index i = l,...,12. To start 
with, let us first construct the local hopping Hamilto- 
nian between the Mn d-orbitals and the p-orbitals of a 
neighboring B ion 'i' at position r^. Let us now take a 
reference frame with the Mn in the origin and the z-axis 
pointing along the direction of this neighboring ion. 
In this reference frame, with a good approximation, only 
the = state of the five Mn d-states hybridizes with 
the Lz = Q state of the B p-orbital. Correspondingly, the 
hybridization between the impurity and this neighbor can 
be approximated as 



V,^V ^tll J] 



(6) 



where d"^ is the annihilation operator for a local Mn 
d-orbital at the origin oriented along the direction n;. 
Similarly, ^'j!. 5. the annihilation operator for the p 
state at the 6 site oriented along the same direction. 
These operators are related to the operators occurring 
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FIG. 2: Fig. a: Local environment of an Mn ion in MgBj. 
The Mn ion is located in a hexagonal cage made of B ions 
represented by large dots. Fig. b: Level structure and crystal 
field splitting of the Mn^""" core states. Two levels (/i — ±2 
and ^ = ±1) are two fold degenerate and the level /i = is 
non-degenerate. The order of the d-levels may depend on the 
details of the cristal field. 
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the -ffo by simple rotations, 



f3=x,y,z 



(7) 

(8) 



where g. refers to states with a quantization axis per- 
pendicular to the B planes, (pi = ~ (pi, and 



{0) = 



2^|sin2, 

'2 V 2 ' 



i (1 + 3 cos 261) 

.1 fl..\^2e 



V 



2 'fsin^l 



(9) 



Summing over all neighboring atoms and expressing all 
operators ^Ti,i3,cr in terms of the band operators, c\^.b,a, 
we then obtain the following hybridization Hamiltonian, 



b,fJ.,(7 



(10) 



where the operator '^b,p,,tT creates an electron with the 
same local d-state symmetry as |/i) in band b, and can 
be expressed as 



6,(7 



U,b (k) 



and V/''^ = VA^^b- In these expressions the normaliza- 
tion factor A^^b has been determined numerically, and is 
defined by the condition that /^ f, (k) be normalized at 
the Fermi surface, 



k 

a* 



^ A 



ie,)e' *'''<efc;a,5. (k) e 



ikRi 



(11) 

(12) 



1 



- J d^k/^,,(k)/;,,,(k) 



(13) 



Symmetry further implies that states belonging to the 
same irreducible representation have the same hybridiza- 

tion: V,^"'' = Vt"'' . 

The above hybridization Hamiltonian generates an ef- 
fective exchange interaction between the Mn spin and the 
conduction electrons in the B bands, since it generates 
charge fluctuations to the Mn^+ and Mn'^^ states. Sec- 
ond order perturbation theory in the hybridization leads 
to the effective exchange Hamiltonian: 



(14) 



6,6',/i,a,/3 



where cr denotes the Pauli matrices, S is the Mn spin, 
and the exchange couplings are given by 



(A')t/(A') 



AE 



A^^bA^^b' I 



(15) 



with AE the characteristic energy of charge fluctuations. 
Note that the symmetry index /x is conserved in Eq.(ll4p. 
thus there are five independent orbital channels of the 
conduction electrons that couple to the impurity spin. 
This is simple to understand on physical grounds: the 
half-filled d-shell has no orbital structure. Therefore, 
a conduction electron that arrives in an orbital state 
must be scattered back to the same orbital channel. How- 
ever, electrons can be scattered between different conduc- 
tion bands, it is only their orbital label that is conserved 
over the scattering process. Therefore, in the absence 
of superconductivity, the channel labels play no special 
role, and the S — 5/2 spin of the Mn ion would be exactly 
screened, resulting in a Fermi liquid statei^ 

By construction, the exchange coupling above satisfy 
J^**' — ^ J^J^^'-, and furthermore, they are equal in 

channels ±/i by symmetry. From Eq. it also follows 
that all the results depend only on a single dimensionless 
coupling proportional to V'^/AE. We define this cou- 
pling as 



(16) 



with Qb the density of states at the Fermi energy in band 
b for one spin direction. Furthermore, in the rest of this 
paper we shall only consider the classical limit, S ^ oo 
with Jj'jj S — finite. In this limit the impurity has no 
dynamics and we can solve the problem exactly. 

To close this subsection let us introduce Nambu 
spinors, $k,6 = {^k'b}'^'^ 



/ Ck,b,T \ 



k,6 



C k,6,i 

„t 



-C 



(17) 



The introduction of these spinors shall simplify our cal- 
culation considerably in the following sections. We can 
rewrite the Hamiltonian in terms of these in a compact 
form. 



H = E<6(ek,b^' + A6r-)*k,6 (18) 

k,6 

+ E lUimlb ^ ■ S $k',6'/,,.6'(k') , 

k.k' ,6,6' 

where the r*'s denote Pauli matrices acting in the pseu- 
dospin (charge) index of the Nambu spinor. In course 
of the derivation we made use of time reversal symme- 
try that imphes /-^,6 (— k) = (—1)^ f^^b i^)^ and doubled 
the Hilbert space so that the components of the Nambu 
spinors in Eq. (|18p must be considered as independent 
variables. 
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III. GREEN'S FUNCTION FORMALISM 



where F (lo) denotes the matrix, 



In this section we shall discuss how the above Hamilto- 
nian can be treated within the Green's function formal- 
ism. In the classical limit the interaction with the im- 
purity in Eq. reduces to a spin-dependent potential 
scattering and, as we show below in detail, the problem 
can be solved exactly. 

In the non-interacting case, J^** — 0, Green's function 
is given by: 



1 



(19) 



and it is a 16 x 16 matrix, diagonal in the band indices. 
In this expression £k,b and A are also diagonal in band 
indices. 

In the presence of impurity scattering we can treat 
the scattering perturbatively, and use multiple scattering 
theory to sum up the series to all orders. The diagram- 
matic expansion of the Green's function is represented 
in Fig. [3l In the first order of perturbation theory the 
self-energy is given by 



s(i) (k, k', c.) = 5] /; (k) I J, (<tS) u (k') 



(20) 



and is independent of the energy to. Here we deliberately 
separated the form factors from the rest of the expression. 
The next order contribution gives 



4''(k,k',^) 



E 



(7^;„(k,q) (^S) 



(21) 



Gl?)(q,c.) (^S) (7,^,;,,(q,k'), 



where we introduced the notation, C^j^, (k, q) = 
1/2/* J (k) J^** /^^f,/ (q). After summing over the momen- 
tum (as explained in Appendix B) and using the orthog- 
onality of the form factors /'s at the Fermi surface we 
end up with the following expression, 



I]Pj(k,k',c.) = 



(22) 



E. 



/;,jk) 



ij^SF(^)ij^S 



/p.b' (k') , 



bb' 



FIG. 3: Diagramatic expansion for the Green's function when 
multiple scattering on the impurity site is considered. The 
solid line represents the full Green's function, the thick line 
represents the non-interacting part of the Green function. 
Each cross represents a scattering on an impurity site, and 
the dotted line stays for the impurity scattering potential, b 
and b' stay for band indices. 



Fbb'ito) 



D 

Sbb'Qb j de 

-D 



Lu — er-' 



(23) 



with D a high-energy cut-oflF that can be removed in the 
end of the calculation. 

Higher order terms can be handled in a similar way. 
The final expression for the Green's function is simply: 



Gbb' (k, k' , c^) = 6kM' Sbb' G];" (k, u) + g[°^ (k, 

E^/;.^« [f,(^) 



/^,,,(k')Gi°)(k',c^). (24) 



By the orthogonality relation Eq. I13[ the quantum num- 
ber fi is conserved. Therefore the T-matrix T(^) can be 
computed independently for each channel /x and is given 
by the following expression: 



f^iuj) = J^S a/2 1 - F(w) J^S a/2 



(25) 



where F{uj) denotes the diagonal matrix F{uj) — F^b'Sbb' ■ 
Note that F {uj) is diagonal in the spin labels. Therefore, 
even order terms in the T-niatrix are spin independent. 
These terms can therefore be referred to as the "charge 
scattering channel" . Odd order terms, on the other hand, 
give spin-dependent contributions and can be referred to 
as a "spin channel" . The even (charge) channel can be 
directly resolved using STM technique while for exper- 
imental observation of the odd (spin) channel contribu- 
tions spin-resolved-STM is needed. 

Impurity bound states and resonances can be iden- 
tified from the pole structure of the T-matrices: True 
bound states correspond to zeros of the determinants 
det{r^^(a;)} on the real axis, and must satisfy \uj\ < At 
for all bands. Zeros in the vicinity of the real axis, on 
the other hand, correspond to resonances. We found that 
each channel generates a bound state, but two of them 
are doubly degenerate by symmetry (/i —p)- 

It is, in general, impossible to find the poles of the 
T matrix analytically, and numerical calculations are 
needed. However, it is generally accepted that, at least 
phenomenologically, superconductivity in MgB2 can be 
explained using a two-band model. With this simple as- 
sumption, the positions of the resonances are given by 
the following equation: 



(1 



g]}a,{±E)){l~gfa^{±E)) 



(26) 



where t;^'' = t^S ^ QbQbJ'/l/' /2 denote the dimensionless 
couplings in channel /i, and 



ab{uj) 



Ab 



1/2 
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In the present case these equations further simphfy due 
to the relation g]}g^^ — {g]^Y to 

gfa,{±E)+gfa2{±E) = l. (27) 

In the hmiting case of J^^ = 0, Eq. (pS)) would give rise 
to two pairs of Shiba states^ for each channel corre- 
sponding to two independent bands. Exchange coupling 
between the two bands, however, removes half of these 
resonances. Similarly, in the realistic situation we thus 
obtain five pairs of Shiba states corresponding to the five 
channels, but two pairs of them are two- fold degenerate 
because of the symmetry j'^ — j'^^. 

IV. DENSITY OF STATES 

Our main purpose is to compute the local tunneling 
density of states (LDOS) and the spin resolved density 
of states near a magnetic impurity for various geometries. 
To obtain a quantitative estimate for the STM spectra we 
performed a lengthy, but straightforward tight-binding 
calculation to determine numerically the form factors 
/^_fc(k), the exchange couplings and the electronic wave 
functions above in various geometries. 

The differential conductivity dl/dV measured by STM 
is proportional with the local density of states which can 
be calculated as the imaginary part of the retarded posi- 
tion dependent local Green's function: 

g,^„(r,c^) = -i-ImTr|G(r,p„,c^)ii^| , (28) 

= ^ E e-*'''-''')i^e^,„,e6.„,G,.,,(k,k',a;). 

Similar to the charge density of states, we can also define 
the spin density of states as 

£»5,a(r,tj) = -^ImTr |G'(r,pa,w) crn ii-^l , (29) 

where n is a unit vector pointing in the direction along 
which we measure the spin density of states. 

The equations above refer to the case where the impu- 
rity is embedded in the bulk. However, both Qs^a and Qs^a 
can be computed easily from the analogue of Eq. ((24)) for 
other boundary conditions too, once the wave functions 
appearing in Eq. ([3|) are known. In the following subsec- 
tions, we first compute the LDOS for an impurity in the 
bulk. Then we study the effect of a semi-infinite half- 
plane with the Mn impurity above and below the first B 
layer. 

A. Impurity in the bulk 

As a first step, we identify the positions of the reso- 
nances for each channel separately from the poles of the 
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FIG. 4: Position of the quasiparticle poles for separate chan- 
nels IJ,. The resonances corresponding to the fi = ±1 chan- 
nels are well resolved for any value of the coupling strength 
g (solid lines). For large enough g, [g > 0.6) the /i = res- 
onance moves inside the gap (dotted lines). The couplings 
corresponding to the /i = ±2 channels are much smaller that 
those in the other channels so the resonances for fi — ±2 
channels are still merged with the superconducting peaks at 
energy A.^ (dashed lines). 

T matrix. In Fig. |4]we show the positions of the bound 
states and resonances obtained as a function of the di- 
mensionless coupling g. The corresponding normalized 
Pz LDOS at the B sites next to the Mn impurity is pre- 
sented in Fig. [5] for different values of g. Due to hexag- 
onal symmetry all B sites around the magnetic impurity 
have the same LDOS. For small values of g the bands are 
slightly interacting and only the most strongly coupled 
fi — ±1 channels give rise to a well resolved resonances 
in the gap for g < 0.4. Increasing the coupling the 
bands are more strongly interacting and the resonances 
corresponding to ^ = channel move inside the gap too. 
This is accompanied, on one hand, by a transfer of weight 
between resonances and secondly by a shift in position of 
each resonance. We also observed small features at en- 
ergies Lu = 7.5meV, i.e. at the energy corresponding to 
Act, due to the coupling between the bands (not shown 
in this figure). 

Fig. [5] also shows the density of states at the next- 
nearest-neighbor sites. The wave functions of the Shiba 
states and thus the amplitudes of the corresponding reso- 
nances in the spectrum depend a lot on the tunneling po- 
sition: The weight and the amplitude of the resonances 
decreases considerably while their position remains un- 
changed. This suppression refiects the local structure of 
Shiba states. At the same time, the coherence peaks near 
the superconducting gap edge gain some spectral weight, 
but they are still quite reduced compared to the bulk. 
Further away from the impurity site the superconducting 
coherence peaks are completely restored and the bound 
states have negligible amplitudes. For generic values of 
the exchange coupling usually two well-separated pairs of 
resonances can be observed, corresponding to the \i = ±1 
and /i = channels. The exchange couplings in channels 
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FIG. 5: Upper panel: Normalized LDOS at the nearest- 
neighbors Boron sites (labeled A in Fig. [2]) for the pz-orbitals, 
for different coupling constants. Lower panel: Normalized 
LDOS at the next-nearest neighbors Boron sites labeled B in 
Fig. [3 



/i = ±2 are much smaller than those in channels /i = ±1 
and /i = 0, and therefore the corresponding bound state 
are merged with the superconducting coherence peak. 




FIG. 6: Spin polarization (odd part of the spectrum and spin- 
up density of states for Site A in Fig. [2] 



The Shiba states are also strongly spin-polarized, as 
is obvious from the spin polarization in the local density 



of states shown in Fig. [6l This fact has an important 
consequence from the point of view of the observability 
of these bound states. As always, a sharp local spectro- 
scopic feature could be difficult to detect if it is over- 
shadowed by the intense continuous background of the 
superconductor. However, the background continuum in 
a superconductor is not, generally speaking, spin polar- 
ized. Thus, even if a Shiba peak happens to be close to 
one of the otherwise dominant BCS coherence peaks, a 
spin-polarized STM can distinguish the Shiba states from 
the continuum^ since the asymmetric part of the spin- 
polarized spectrum has sharp peaks at the resonances 
but is predicted to be featureless otherwise. Therefore 
spin-polarized STM is clearly an ideal tool to identify 
the multiple Shiba states. 



B. Impurity in the vicinity of a surface 

As we mentioned already, the effect of a surface can 
be taken into account by simply modifying the wave 
functions that appear in the expansion of the operators 



(30) 



Here kj^ is the in-plane momentum and is the momen- 
tum perpendicular to the surface. Note that the surface 
breaks translational symmetry along the z direction, and 
therefore only fcz > values are permitted. The wave 
functions above must satisfy the appropriate boundary 
conditions, and can be expressed within our tight bind- 
ing formalism as 



Jk±R_ 



-V2 sin {k^Z) , 



(31) 

with Z = corresponding to the first layer in the vacuum. 

Our calculatuions for an impurity in the bulk can easily 
be extended to this case as well with minor modifications. 
If the magnetic impurity is well inside the bulk, we re- 
cover the results discussed in the previous subsection. In 
Fig. [7] (upper panel) we represent the LDOS at nearest 
neighbor B atoms for the case when the Mn impurity is 
below the top B layer. The amplitudes of the resonances 
are slightly reduced in this case compared to the bulk sys- 
tem and also the positions are modified due to the local 
density of states that is slightly modified in the vicin- 
ity of the surface. Moving away from the impurity, the 
weights of the resonances start to decrease and the su- 
perconducting coherence peaks are gradually recovered. 
In this configuration, at sites more than two lattice con- 
stants away from the impurity site the superconducting 
coherence peak is already completely recovered. 

For spin-resolved scanning tunneling spectroscopy, the 
tunneling current can be separated into an unpolar- 
ized part /q, which depends only on the LDOS, and a 
spin-polarized contribution Ip given by the projection 
of the local magnetization density at the tunneling site 
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FIG. 7: Upper panel: Normalized density of states for site 
A as function of frequency for different values of g in the 
geometry when the Mn impurity is below the first B layer. 
Lower panel: Normalized LDOS at the next-nearest neighbors 
Boron sites B of Fig. O when the Mn impurity is just below 
the first layer 



onto the magnetization direction of the tip. The spin- 
polarized contribution to the local differential conductiv- 
ity is therefore proportional to the magnetization density, 
dIp/dV oc Pt cos9gs{ri,uj = eV), where Pt denotes the 
polarization of the tip, is the angle between the mag- 
netization axes of the tip and the impurity spin. 

In Fig. [8] we present the the local spin polarization at 
site A for g = 0.491. For the same reasons as before, 
only the contribution of the Pz orbitals is shown. The 
relative orientation of the impurity spin and the tip can 
also be fixed by a small external magnetic field in these 
experiments. However, the angle 9 is not arbitrary even 
in the absence of an external field, since in the vicinity 
of a ferromagnetic STM tip a magnetic impurity would 
be presumably aligned with the magnetization of the tip 
due to stray fields. The most important feature we ob- 
serve is a transfer of weight from states in the gap to 
states in the continuum due to the inter-band coupling 
through the magnetic impurity. The inset presents the 
total spin-polarized tunneling density of states for the 
same coupling (7, for a complete polarization of the tip, 
Pt = 1 and a perfect alignment, 9 = 0. 



n=+/-i 



g=0.491 



^=o 



-5.0 -2.5 0.0 2.5 5.0 



0.0 
fo(meV) 



FIG. 8: Spin polarization and spin resolved density of states 
for Site A of Fig. [2] 



V. CONCLUSIONS 

We presented a detailed theoretical investigation of the 
effect of a single Mn magnetic impurity on the supercon- 
ducting properties of MgB2. Our description is based 
on a microscopic model which assumes nearest neigh- 
bors hopping from the localized orbital of the Mn to the 
neighboring B orbitals. We have shown that a magnetic 
impurity generally induces multiple Shiba states in the 
electronic structure of MgBj. In particular, for Mn we 
found five pairs of Shiba states in the gap, two of which 
were two-fold degenerate. We have taken into account re- 
alistic band structure and the efi^ect of surface states on 
the local spectrum. Our calculation of both conventional 
and spin-resolved STM^^ spectra near the impurity site 
showed that these states can be clearly resolved by both 
methods. Similar multiple Shiba states should appear 
in other superconductors due to the internal structure of 
the magnetic impurity. 

It is intriguing to speculate what these local probes will 
eventually see in an actual experiment. Clearly, despite 
decades of pioneering investigation, local spectroscopy of 
spin impurity states in a superconductor still has the po- 
tential of revealing new features that have not yet been 
documented. For example, our calculations assume clas- 
sical spin degrees of freedom, whereas the experimental 
measurements could reveal - besides a classical behavior 
- effects of screening of a quantum spin by the supercon- 
ductor, leading to either a full screening or a reduction 
of the effective spin carried by the impurity. The quan- 
titative discussion of such effects goes beyond the scope 
of the present paper. Nevertheless, our calculations will 
provide an important benchmark for comparison with ex- 
periments, a benchmark that includes, for the first time, 
the presence of multiple channels of scattering. 

The results we obtained in this paper are relevant 
and relatively easy to generalize for other compounds. 
For example, recent STM measurements have focused 
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on Ti impurities in another multiband superconductor 
Sr2Ru04. While these experimental results are prelim- 
inary as the magnetization state of Ti is not clear, and 
there are several differences between MgB2 and Sr2Ru04, 
it is clear that our framework provides a suitable plat- 
form for studying Sr2Ru04 as well. As we mentioned 
above, MgB2 crystallizes in the hexagonal AlB2-type 
structurefiS, and the band structure of MgB2 is also some- 
what peculiar. Nevertheless, as shown in the seminal pa- 
per of Nozieres and Blandin-, although the form of the 
exchange Hamiltonian depends a lot on the specific ma- 
terial and point group considered, in most cases, similar 
to Mn-doped MgB2, several channels of conduction elec- 
trons couple to the local impurity degrees of freedom, 
and result in multiple Shiba states. Therefore, that the 
appearance of multiple Shiba states is a rather general 
phenomenon. 

An interesting result of our analysis is that, although 
it may be difficult to resolve a Shiba state close to the 
coherence peak with conventional STM methods, the an- 
tisymmetrical part of a spin-resolved STM clearly sepa- 
rates these states in the STM spectrum. The weight of a 
given pair of Shiba states may, however, be very sensitive 
to the particular atomic state into which electrons tun- 
nel from the STM tip, and depends also on the precise 



-ffa;,l;a;,l(k) 
-ffy,l:y,2(k) 



All the other matrix components of the Hamiltonian ma- 
trix can be written in terms of those given in Eq. (|A1() 

as follows: Hy^i-y^iCk) ^ 7?x,2:x,2(k) = iJy,2;y,2(k) = 
Hx,l;x,l{'^), ^^a,l;K,2(k) — HxS;y,20<-), -ffz,2;z,2(k) ~ 

Hz.i-z.iO^)- AH the other elements are equal to zero. This 
matrix is Hermitian, Hji = H*^ for i ^ j- The best fit to 
other calculated band structurei^ is obtained for the fol- 
lowing set of parameters: exy = — 8.6eV, = — 1.5eV, 
t = 2.0 eV, tz = 2.5 eV, t\\ = 4.5 eV, t±_ = 1.8 eV, 
txy = 0.1 eV. In our calculation, t and are the hop- 
ping integral corresponding to the pz orbitals: t is the 
in-plane hopping between the nearest neighbors (tt bond- 



position of the tip. 
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APPENDIX A: TIGHT-BINDING 
HAMILTONIAN 

In this Appendix we present the basic results obtained 
from the tight-binding analysis of the bulk system. The 
matrix elements of the Hamiltonian given in Eq. (|4]) are 
given by: 



(Al) 



I 

ing), and tz is the out of plane hopping (ct bonding). 
The parameters i|| and tj_ denote a and 7r-like hopping 
integrals for the in-plane Px,y orbitals. Finally, the out 
of plane hopping integral is given by txy This param- 
eter is very small, so there is practically no dispersion 
along the T — A line. The corresponding density of states 
(DOS) has been calculated in the framework of Green's 
function formalism as g{uj) = X^k^™^'-''' (^''^) where 
Gb (k, Lu) is the Green's function corresponding to every 
band b. The resulting DOS is presented in Fig. [9l The 
values for the DOS at the Fermi surface for the bands 
that cross the Fermi surface are Qx — 0.081 states/eV, 



€xy '^ixy cos 



t± + -t» + -t± ) exp ( —i—ky I 2 cos 



V3 



t|| — <j^) exp ^i^iky I 2i 



I sm - 



V3, 



h + \-t 



1 



.4 4 

€z — 2tz COS 



t|| exp -i-ky 2cos—kx , 



/ /V3, 3, 
t 1 + exp —I -r-fcx + 7:^ 
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In the numerical calculations we used a 100 x 100 x 100 
discretization of the first Brillouin zone. For each site 
in the discretized lattice we calculated the energy values 
corresponding to band b. We then tested if one of these 
10^ cells overlapped with the shell of thickness de. If it 
did, we generated a mesh of 15 x 15 x 15 within this cell 
to compute the cell's contribution to Eq. IB4I 

In our calculations, the number of points around the 
Fermi surface was larger than 10^ within an energy shell 
of lOmeV. This was used to evaluate the average of the 
form factors and their Fourier transforms at the Fermi 
surface with a precision of ^ 10^"^. 



APPENDIX C: MOMENTUM SUMMATION 



FIG. 9: Top: Density of states for px,y and Pz bands. Bot- 
tom: Total density of states. The vertical line represents the 
position of the Fermi energy. 

Qy = 0.13states/eV, g^i = gz2 = 0.75 states/eV, in rea- 
sonable agreement with more sophisticated band struc- 
ture calculationsii. 



In this section we explain the method that we used 
to evaluate the momentum summation in the first Bril- 
louin zone. During the calculations we have to evaluate 
expression of the form: 

l^<^(k)Gi°)(k,c.) , (CI) 



APPENDIX B: AVERAGE OVER THE FERMI 
SURFACE 

Throughout our analysis we have to evaluate averages 
over the Fermi surface. For a given band we have to 
calculate 



(Bl) 



where Sb represents the Fermi surface area for band b and 
if (k) is a momentum dependent function. The Fermi sur- 
face was obtained in our calculation by numerically solv- 
ing the equation ek,6 = 0. To evaluate (|B1[) we replace 
the integration over the Fermi surface with an integra- 
tion over an energy shell of thickness de. First the area 
of the Fermi surface can be calculated as 

Sb^ f d^k=^ f d^klVek.bl . (B2) 



de 



where ip (k) is a momentum dependent function (usu- 
ally the form factor or a combination including form 
factors and other momentum-dependent functions) and 

G^"^ (k, uj) is the free Green's function. The free Green's 
function depends on momentum only through the energy 
of the given band £k,6- We approximated therefore the 
summation as 



i5:^(k)Gr (k. 



(C2) 



oo 

Qb J deGf\e,u;)^J^ d^]^^{\^) , (C3) 



with Qb the density at the Fermi surface in band b. The 
integration of the Green's function over the energy can 
be done analytically and the result is simply 



hell 



In a similar way the average of any momentum- 
dependent function can be evaluated as: 

f ^{k)d^k=^ f d^k\Weu,b\v{k) . (B3) 

J St J shell 

Our quantity is therefore given by the expression: 



1 / ,.rkw^k./--rf^k|vsk,|^(k) 

Js, 



Isheii rf'k|Vek,&' 



Fbitu) 



Qb 



de 



uj — eT'^ 



AbT"" 



(C4) 



For |a;| < A{, the function Fb{uj) has only real parts and 
simplifies to Fb{u;) = -iTQb{w + AhT^)/(A2 - tj2)-i/2^ 
while for > Af, it is purely imaginary. The other 
term which represents an average over the Fermi surface 
was calculated numerically as explained in Appendix B. 
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